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i INTRODUCTION 


It  is  clear  that  a fundamental  understanding  of  turbulence  would  be 

of  great  practical  significance.  Since  it  is  suspected  that  on  the  very 

smallest  scales  turbulence  has  certain  universal  features,  it  might  be 

hoped  that  these  would  be  revealed  by  a direct  numerical  simulation. 

Unfortunately  the  problem  of  fully  developed  turbulence  is  one  involving 

very  nu  .y  degrees  of  freedom.  Indeed,  if  R is  the  Reynolds  number , the 

necessary  degrees  of  freedom  to  be  treated  in  a given  dimension  are  of 
3 '4 

the  order  of  R .In  three  dimensions  this  number  is  then  of  the  order 
9 /4 

of  R (see,  fo»  example,  Landau  and  Lifshitz,  1959:  for  more  detailed 

* 

estimates,  see  Case  et  al . , 1973).  Since  time  steps  are  restricted  by 

3 

space  steps,  the  number  of  needed  calculations  then  grows  as  R (Case 
et  al.,  1973).  For  the  large  values  of  R occurring  in  practice  the 
computations  then  become  prohibitive  with  present  (or  soo. ’-to-be-available) 
computers . 

Under  these  circumstances  it  is  perhaps  reasonable  to  back  off  a bit 
and  ask  what  can  be  done  in  the  way  of  a numerical  investigation  of  the 
transition  from  laminar  to  turbulent  flow,  here,  unless  otherwise  1 ndi 
cated,  we  restrict  ourselves  for  simplicity  to  the  incompressible  case — 
i.e.,  our  flows  are  described  by  the  Navier-Stokes  equations 

~ + v . Vv  = -VP  + vV  v (1) 

at  ~ 

V • v = 0 . (2) 

* 

References  are  listed  at  the  end  of  the  report. 
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Knowledge  of  such  transitions  also  has  considerable  significance  for 
practical  designs.  The  hope  is  that  we  will  be  able  to  treat  such  prob- 
lems with  presently  available  computers,  since: 

(1)  Transition  frequently  takes  place  at  relatively  low 
Reynolds  number. 

(2)  Transitions  are  frequently  seen  to  develop  through  a 
number  of  stages  in  which  the  flow  is  relatively  large- 
scale  (i.e.,  involving  relatively  fev  degrees  of  freedom) 

Experimentally,  the  transition  takes  place  in  many  different  ways. 
Indeed,  for  a given  practical  situation  more  than  one  mechanism  may  be 
responsible.  There  are  a number  of  ways  in  which  one  can  characterize 
the  different  possibilities.  None  of  these  is  completely  satisfactory 
in  that  the  descriptions  are  not  necessarily  unique  or  mutually  exclusive. 
One  such  cataloging  is  in  terms  of  whether  the  transition  occurs  in  a 
flow  that  is  inviscidly  stable  or  not.  Another  rather  convenient  dis- 
tinction is  between  transitions  in  "free'’  or  "bound"  boundary  layers 
(Sato  and  Kuriki , 1961).  in  essence,  the  difference  is  between  flows 
that  become  turbulent  rapidly  and  those  that  become  turbulent  gradually. 

It  is  the  latter  class  that  we  expect  to  be  most  amenable  to  numerical 
simulation.  However,  as  we  will  see,  the  division  is  by  no  means  sharp. 

Accordingly,  we  consider  .he  calculation  of  transition  in  a variety 
of  simple  situations  that  have  shown  different  forms  of  transition.  We 
start  with  problems  in  which  there  is  much  detailed  experimental  infor- 
mation available  and  some  theoretical  insight.  As  will  be  seen,  this 
knowledge  is  very  useful  in  showing  us  what  size  calculation  must  be 
performed  and  what  one  expects  to  be  able  to  describe.  Assuming  that 
the  initial  calculations  are  successful,  we  then  indicate  others  that  can 
be  done  that  hopefully  will  lead  to  methods,  in  which  there  is  confidence, 
for  calculating  transition  in  situations  of  practical  interest. 
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Our  program  is  then  as  follows.  First  we  consider  various  situations 
in  which  transition  could  be  calculated.  These  are  enumerated  essentially 
in  the  order  of  what  we  consider  their  priority.  m each  case  we  describe 
the  physics  of  the  situation  and  our  present  understanding.  In  various 
degrees  of  detail  we  discuss  the  limitations  of  the  calculations  (reso- 
lution needed  and  uncertainties).  An  estimate  is  then  made  of  the  compu- 
tation effort  required.  Second,  we  give  some  recommendations  as  to  how 
a program  to  do  these  things  could  be  implemented.  These  include  cost 
per  year,  number  of  years,  organization,  and  supervision. 

It  is  felt  certain  that  such  a program  would  contribute  significantly 
to  the  present  state  of  the  art  with  a relatively  modest  expenditure  of 

funds.  It  is  hop~d  this  will  lead  to  a capability  of  doing  some  practi- 
cal calculations. 
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II  THE  FLAT  PLATE 


A 


* 


> 


As  a first  calculation  we  suggest  the  numerical  study  of  the  transi- 
tion in  the  boundary  layer  in  the  flow  over  a flat  plate.  The  advantages 
of  this  are: 

fl)  Codes  ior  this  have  been  developed  and  some  calculations 
have  been  performed  (Grosch,  1974;  Orszag,  1974). 

(2)  Excellent  detailed  experimental  information  is  available 
(Klebanoff  et  al.,  1961;  Kovasznay  et  al,,  1962), 

(3)  Considerable  theoretical  insight  is  available.  Thus 
the  initial  stage  of  laminar  instability  is  well  under- 
stood (Schlichting , 1960).  Theoretical  models  give  a 
qualitative  picture  of  what  happens  in  the  later  stages 
(e.g.,  Benny,  1961;  Stuart,  1962;  Lin  and  Benny,  1962; 

Greenspan  and  Benny,  1963). 

There  are  also  disadvantages: 

(1)  As  indicated  in  our  earlier  report  (cf.  Schlichting, 

1960,  p.  386)  there  are  some  ambiguities  concerning 
downstream  boundary  conditions. 

(2)  At  the  last  stages  before  complete  turbulence  is 
obtained,  small-scale  motions  are  found.  There  is 
some  question  as  to  whether  a calculation  designed  to 
follow  the  earlier  stages  of  transition  can  resolve 
these . 

We  return  to  these  two  points  later. 

To  illustrate  what  is  involved  in  the  calculation  of  transition  of 
the  boundary-layer  flow  over  a flat  plate,  we  give  a brief  (and  idealized) 
version  of  the  experiment  described  by  Klebanoff  et  al.  (1961). 

Over  a flat  plate  located  at  x 2 0,  y = 0 there  flows  a fluid  with 
velocity  Uoo  in  the  x direction  toward  the  plate.  The  pressure  gradient 

Preceding  page  blank  , 


a 


ma 


I 


i 


is  adjusted  to  be  zero.  A ribbon  is  located  above  the  plate  in  a span- 
wise  (z)  direction  perpendicular  to  the  flow  at  some  distance  downstream 
from  the  leading  edge.  Under  the  ribbon  there  are  strips  uniform  in 
length  and  uniformly  separated.  With  the  ribbon  stationary  the  x component 
(u)  of  the  velocity  is,  for  any  x and  z,  as  shown  in  Figure  1. 
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FIGURE  1 MEASURED  VELOCITY  DISTRIBUTION  IN  THE  LAMINAR  BOUNDARY  LAYER 
ON  A FLAT  PLATE  AT  ZERO  INCIDENCE.  Source:  Schlichting  (1960). 

The  ribbon  is  then  vibrated  at  a frequency  f.  At  small  amplitudes 
of  oscillation  the  results  are  as  follows.  For  a fixed  distance  x down- 
stream the  velocity  varies  essentially  sinusoidally  in  time  and  with  wave 
length,  which  is  the  periodicity  of  the  strips  in  the  z direction.  The 
average  velocity  distribution  as  a function  of  y retains  essentially  the 
form  of  Figure  1.  As  we  go  downstream  the  amplitude  of  the  oscillations 
in  u first  increases  and  then  gradually  decreases. 
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At  , larger  amplitude  ef  oscillation  the  picture  changes.  As  „e 
proceed  downstream  from  the  ribbon  the  situation  Is  first  os  above. 

Further  down,  however,  there  develops  a pronounced  three-dimensional 
Structure  characterised  by  large  spanwlse  variations  in  wave  amplitude 
«th  peaks'-  and  "valleys  " occupy  ng  fixed  spanwlse  positions.  Associated 
•ith  this  variation  In  amplitude  there  is  also  a spanwlse  variation  In 
local  mean  velocity  such  that  there  is  a defect  at  a peak  and  an  excess 
at  a valley.  Further  downstream  there  occurs  an  abrupt  increase  in 
amplitude  at  a peak  that  is  characterized  by  a series  of  Intense  low- 
velocity  pulses  evidenced  by  "spikes"  In  an  oscillogram  of  strea.wise 
fluctuating  velocity.  At  first  a single  spike  appears  for  each  cycle  of 
the  primary  oscillation.  The  spikes  Increase  in  number  as  we  go  down- 
stream and  ultimately  blend  into  fully  developed  turbulence. 


Present  theoretical  understanding  is  the  following.  In  the  absence 
of  the  vibrating  ribbon  the  stationary  laminar  flow  is  adequately 

described  by  the  boundary- layer  equations  for  the  velocity  [u(x ,y) ,v(x ,y) , 
w = 0]: 


udu  vdu  _ $ u 

3x  3y  - V 2 

dy 


(3a) 


3u  6v 

r~  + — = o 

OX  0y 


(3b) 


u-v=0aty=0  , limu=U 


y-+oo 


Those  equations  have  the  Blasius  similarity  solution 


u = U f (T| ) 


v = 


f (71) 

i x i' 


(4) 


where 


Tl  = 


With 


f (Tj)  = cp 7 (Tl)  , r 1 (Tl ) = 1/2 [ Hep / (T| ) - cpl 


we  have 


epep ,/  + 2cp  = 0 


(5) 


subject  to 

cp(O)  = cp'(O)  = 0 , cp'(»)  = 1 

Numerical  integration  then  gives  the  velocity  profile  of  Figure  1 — in 
excellent  agreement  with  experiment. 


Remarks : 

(1)  Two  conventional  measures  of  the  boundary- layer 
thickness  used  are: 

(a)  6 = that  value  of  y for  which  u = 0.991^, 

Numerically, 


<b)  M (1-t)d’'*‘1-73^  ■ 

0 

(2)  For  the  experiment  of  Klebanoff  et  al.  (1961)  approxi- 
mate numbers  are  U®  ~ 1500  cm/s,  v = 0.15  cm^/s,  and 

x ~ 100  cm.  Then,  5 j ^ - 0.5  cm.  Typically 

ribbon  ribbon 


8 


the  variations  described  above  occurred  within  40  cm 
of  the  ribbon.  At  this  point  6 has  changed  by  only 
~20%. 

(3)  From  Figure  1 and  the  above  we  conclude: 

(a)  The  velocity  distribution  varies  smoothly  within 
the  boundary  layer. 

(b)  After  but  a few  6 from  the  plate  (~1  cm  for  the 
experiment)  the  flow  is  essentially  arbitrarily 
close  to  the  free-stream  velocity. 

Tc  study  the  effect  of  external  disturbances  it  is  natural  to  turn 
to  linear  stability  theory  Thus,  in  Eqs  . (1)  and  (2)  we  substitute 

u = u + 6u  (7) 

~ ~o  ~ 

where  u^  is  the  original  laminar  flow,  and  we  retain  only  terms  linear 
in  the  perturbation  6u.  Assuming  that  a time  dependence  6u  ~ e i<Ut  leads 
to  an  eigenvalue  problem.  If  it  is  found  that  Im  m = ^ is  >0,  then 
initial  perturbations  can  grow,  and  we  have  instability.  In  practice, 
another  approximation  is  made.  Stability  locally  is  studied  by  assuming 
that  the  profile  (y  dependence)  is  such  that  the  one  at  a given  position 
*o  really  extends  from  -<»  < x < ».  (Since  in  the  problems  considered 
here  the  true  x variation  of  the  profile  is  very  slow,  this  is  a readily 
justifiable  approximation).  In  this  case  th?  eigenvalue  problem  is 
simplified.  Since  there  is  no  explicit  x-dependence  we  can  limit  our- 
selves to  investigating  solutions  with  x and  z dependence  of  the  form 
-(orx+Pz) 

~ e . [Since  Squire's  theorem  (Schlichting , 1960,  p.  386)  assures 

us  that  the  onset  of  instability  occurs  first  when  p = 0,  the  original 
solutions  of  the  eigenvalue  problem  were  obtained  for  this  case].  With 
this  form  of  x-dependence  the  eigenvalue  problem  can  be  regarded  as  one 
in  which  or  is  given  (real)  and  the  complex  uj  is  sought. 

Figure  2 gives  the  results  for  an  old  calculation  of  the  curve  for 

* 

neutral  stability  in  the  or6  , R = U 6*/v  plane.  It  may  be  noted  that 

0 * 09 
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FIGURE  2 CURVE  OF  NEUTRAL  STABILITY  FOR  THE 
WAVELENGTHS  a8*  OF  THE  DISTURBANCES 
IN  TERMS  OF  THE  REYNOLDS  NUMBER  R 
FOR  THE  BOUNDARY  LAYER  ON  A FLAT 
PLATE  AT  ZERO  INCIDENCE  (Blasius  Profile). 
Source:  Schl ichting  (1960). 


for  instability,  R must  be  greater  than  a critical  value  R ~ 500 

* ° **C 

and  06  must  be  less  than  about  0.36  (i.e.,  there  is  a minimum  wavelength 
* * 

min  ~ 2"6  /0‘36  ~ 17’56  ~ 06  for  instability).  For  the  experiment 

described  above  this  wavelength  is  X ~ 3 cm,  while  the  distance  x 

min  c 

downstream  from  the  leading  edge  where  instability  can  possibly  occur  is 

x ~ 80  cm. 
c 

For  the  experiment,  however,  it  is  somewhat  more  relevant  to  con- 
sider a real  uj  and  find  the,  in  general,  complex  o.  (With  a = Im  a < 0 , 
we  have  spatial  amplification.)  A neutral  curve  in  the  u>,  R plane  is 
shown  in  Figure  3.  We  see  that  there  is  a maximum  frequency  for  which 
such  spatial  amplification  holds.  This  is 


* * ' ’ **  ** 

!U 3* 

FIGURE  3 CURVE  OF  NEUTRAL  STABILITY  FOR  THE 
FREQUENCY  /?r  OF  THE  DISTURBANCES 
AND  WAVE  PROPAGATION  VELOCITY  cr 
FOR  THE  BOUNDARY  LAYER  ON  A FLAT 
PLATE  AT  ZERO  INCIDENCE  (Blasius  Profile) 
Source:  Schlichting  (I960). 


For  the  experiment  described,  this  is; 


f ~ 300  c/s 
max 


We  note  also  from  Figure  3 that  the  -Hoc 

b hat  the  phase  velocity  (c  = m/a)  of  the 

neutral  waves  is  of  order  0.4U  or  less 

00  * 


TTe  theoretical  picture  of  the  early  stages  of  the  experiment 
“escribed  is  theu  the  folio**.  The  oscillating  ribbon  excites  eaves 
of  wavelength  X > 3 on.  These  are  then  amplified  as  we  go  downstream. 

the  initial  perturbation  is  sufficiently  small  we  will  pass  through 
the  right-hand  neutral  curve  before  the  altitudes  are  large  enough  to 
near  assumption.  As  we  go  farther  downstream  the  disturbanc 
will  then  be  damped  out.  On  the  other  hand,  if  the  initial  disturbance 

e"°Ugh  “ WU1  ^ ■*“«-  «tM.  the  unstable  region  down- 

stream  to  a point  where  linear-stahi  n tv  + », 

near  stability  theory  no  longer  applies  and 

cannot  be  used  for  prediction 


h 

i t 
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One  further  result  ef  linear-stability  theory  of  importance  for  us 
“ the  Sh“PC  o£  th<!  eirenf unct Ions  corresponding  to  the  unstable  mode. 
From  Figure  4 we  see  that  they  vary  rapidly  in  the  part  of  the  boundary 
layer  between  y - 0 and  , . 0.2  - 0.46.  (For  the  experiment,  this  number 
is  y - 0.1  - 0.2  cm.)  Outside  this  region  the  variation  is  rather  smooth. 


FIGURE  4 VARIATION  OF  AMPLITUDE  OF  THE  u'-FLUCTUATION  FOR 
TWO  NEUTRAL  DISTURBANCES  IN  A LAMINAR  BOUNDARY 
LAYER  ON  FLAT  PLATE  AT  ZERO  INCIDENCE.  The  curves 
labeled  I and  II  correspond  to  the  two  neutral  disturbances  I and 
■*  m Fi9ure  2.  Source:  Schlichting  (1960). 

The  further  development  of  the  flow  toward  turbulence  is  somewhat 
less  know,,  theoretically.  From  models  (Benny,  1961;  Stuart,  1962;  Lin 
and  Benny,  1962;  Gre  nspun,  .1962;  Greenspan  and  Benny,  1963)  and  experi- 
ment (Klebanoff  et  al.,  1961)  it  appears  that  the  amplified  two-dimensional 
Tollmien-Schli chting  wave  interacts,  as  a result  of  nonlinear  terms, 
with  true  3-dimensional  disturbances  produced  by  the  spacers.  This 
combination  then  gives  rise  to  the  distorted  boundary  layer  which  in 
turn  results  in  the  production  of  the  high-frequency  disturbances.  While 


1 
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the  primary  contribution  of  a numerical  simulation  would  be  to  allow  the 
development  of  the  nonlinear  interaction,  two  qualitative  conclusions 
can  be  drawn: 

(1)  Resolution  in  the  spanwise  (z)  direction  must  be 
sufficient  so  that  the  expected  peaxs  and  valleys 
can  be  adequately  differentiated. 

(2)  While  the  Tol lmien-Schlichting  waves  have  their 
largest  amplitudes  and  sharpest  variation  close  to 
the  plate  (~0.l6)  the  fluctuations  near  breakdown 
have  their  largest  values  much  farther  out  (~0.56) — 
cf . Figure  5.  From  this  we  conclude  that  even  at 
these  distances  the  resolution  cannot  be  too  coarse. 


FIGURE  5 DISTRIBUTION  OF  INTENSITY  OF  u-FLUCTUATION  ACROSS  BOUNDARY 
LAYER:  145-c/s  WAVE,  U ,/v  = 3.1  X 105  fH.  Source:  Klebanoff  et  al. 

(1961). 


Let  us  now  turn  to  the  numerical  simulation  of  the  above  experiment. 
Clearly,  finite  computer  capacity  restricts  us  to  a discussion  of  the 
flaw  to  some  finite  region  of  space  and  time.  In  the  x (streamwise) 
direction  this  should  begin  somewhat  downstream  from  the  leading  edge  of 
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the  Plate  to  assure  the  Initial  establishment  of  the  laminar  flow  but 
upstream  from  the  ribbon  to  diminish  spurious  effeets  due  to  Ineorreet 
boundary  conditions.  A reasonable  guess  .ould  be  to  star,  at  a position 
‘2-uM^n  fl-°”  “e  ribbon  (where  X^  is  the  minimum  unstable 

lollmien-Schlichti ng  wavelength).  Downstream  from  the  ribbon  we  would 

ncludt  all  positions  till  the  onset  of  turbulence.  Experimentally 
this  is  -SX^.  However,  as  indicated,  later  ambiguities  the  boundary 
conditions  suggest  that  the  end  position  be  chosen  as  far  downstream  as 
feasible.  A compromise  might  be  that  the  region  In  a to  be  described 
is  o,  the  order  of  (10-20, X^.  Since  in  the  span.ise  <„  direction  we 
are  describing  a periodic  behavior,  we  can  restrict  ourselves  to  a dis- 
tance equal  to  one  wavelength  of  the  spa., wise  perturbation.  In  the 
direction  perpendicular  to  the  plate  we  can.  since  the  laminar  flow  and 
Its  perturbations  go  to  zero  rapidly  outside  the  boundary-layer  thlchness 
, restrict  ourselves  to  a distance  of  a few  6 (say,  26).  tl«  we 

would  like  to  follow  the  flow  for  some  significant  number  of  oscillations. 

If  T be  the  period  of  the  ribbon  the  time  of  interest  may  be  of  the  order 
of  (5-10)T. 

NOW  what  do  these  dimensions  imply  about  the  number  of  calculations 

to  be  performed?  The  problem  is  to  solve  Eqs.  (1,  and  (2,  in  the  region 

described,  subject  to  suitable  boundary  conditions  (discussed  below,. 

To  estimate,  we  imagine  the  calculation  done  by  finite  Hi  f r 

1 Dy  unite-difference  methods 

ln  =U  dlm0nS1°“-  (If  ~ ■»  -ed.  a rule  of  thumb  might 

he  that  the  number  of  modes  in  , given  direction  be  -(1/2  - 1/3,  of  ,he 

number  of  grid  points.  the  -direction  we  want  to  describe  a distance 


Of  the  order  15X  . Assuming  that  this  Is  adequately  resolved  by  8 


min 


points  per  wavelength  we  have  something  like  160  points  in  the  x-dlrectlon. 
i.o.,  in  the  notation  below.  L ~ 160.  In  the  z-dlrectlon  we  have  a 
distance  of  only  one  perturbation  wavelength,  but  to  have  a chance  of 
describing  the  variation  near  breakdown  we  need  a fine  resolution.  A 
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guess  is  16  points  in  this  direction  (N  = 16).  In  the  6-direction  we 
have  a total  distance  26,  but  as  we  have  seen  in  the  region  0 £ y < 0.26, 
there  are  very  rapid  variations.  Hence,  with  finite— difference  method® 
one  would  use  either  some  coordinate  mapping  or  a varying  grid  size. 
Assuming  the  latter,  one  might  use  of  the  order  of  20  points  for  the 
first  0.36  and  20  irore  for  the  remainder — i.e.,  of  the  order  of  40  grid 
points.  The  sizes  of  the  time  steps  are  essentially  determined  by 
stability  requirements — i.e., 

vAt  < Ar 

where  v is  some  appropriate  velocity  and  Ar  is  some  mesh  spacing.  Since 

U is  the  largest  velocity  and  it  is  in  the  x-direction,  a reasonable 
00 

criterion  (Grosch,  1974)  might  be 


. „ Ax 

Lt<~ 


With  the  numbers  suggested  above,  this  is 


0.375  -4 

At  < : ~ 2.5  X 10 

1.5  X 10 


-3 

For  T = 1/150  c/s  = 6.7  X 10  , this  implies  of  the  order  of  30  time 
steps  per  cycle  and  therefore  a number  of  time  steps  Nt  ~ 300  to  500. 

Computer  requirements  based  on  these  numbers  are  discussed  below, 
First,  however,  we  include  a few  caveats. 


A.  Boundary  Conditions 

A well  posed  problem  for  solutions  of  Eqs . (1)  and  (2)  would  be 
one  in  which  the  initial  velocity  is  prescribed  everywhere  within  the 
region  described  and  the  velocity  is  given  everywhere  on  the  boundary 


r 


r 

r 


4 
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for  all  time  (subject  to  J n • u ds  = 0) . Unfortunately  we  do  not  know 

this  information.  At  y = 0 we  do.  There  all  components  of  the  velocity 

« 

are  zeio.  At  y = 26  it  is  clearly  correct  to  a very  high  degree  of  ap- 
proximation to  take  the  velocity  as  that  of  the  free  stream.  At  the 
upstream  edge  of  our  region  it  is  reasonable  to  assume  that  the  velocity 
is  of  the  Blasius  form  plus  whatever  turbation  one  might  want  to  put 
in.  In  the  spanwise  direction  a periodicity  requirement  is  reasonable. 
What  we  do  not  know  is  the  velocity  at  tic  downstream  edge  of  our  region. 
Two  somewhat  ad  hoc  suggestions  have  been  made.  Grosch  (1974)  suggests 
"extrapolating”  from  the  velocity  inside.  If  L is  the  downstream  edge 
of  the  region,  this  formally  is  of  the  form 

fL 

u(L,y,z,t)  = I u(x,y,z,t)  g(x)  dx  + 

J 0 

In  some  sense  this  may  be  regarded  as  a kind  of  radiation  condition, 
which  might  indeed  be  reasonable.  How  well  the  actual  choice  made  of 
g (x)  and  h (x)  approximates  such  a condition  is  not  clear.  Also,  there 
is  an  uncomfortable  situation  in  that  the  pressure  at  x = L is  calculated 
from  du/dt , while  to  advance  forward  in  time  du/dt  is  calculated  using 
the  resulting  p.  The  possibility  of  an  instability  occurs. 

Orszag  (1974)  proceeds  by  dropping  certain  terms  as  small.  With  the 
emasculated  equations  he  demonstrates  uniqueness  if  at  all  boundary 
points  either  p or  v • n is  given  as  well  as  v for  points  such  that 

~ * "out  < °*  (Noto:  this  demonstration  involves  an  additional,  but 

reasonable,  approximation.)  The  difficulties  here  are: 

(1)  We  do  not  know  p or  v . n on  the  downstream  boundary. 

(2)  We  do  not  know  v so  that  we  cannot  even  determine  when 

v . n <r  o. 

~ -out 

(3)  The  approximations  involved  need  to  be  justified.  (In 
fairness  it  must  be  noted  that  Orszag  has  recently 
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du 

7”;  (x,y,z,t)  h (x ) dx 
01 
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suggested  sc  me  modified  boundary  conditions  about  which 
we  have  insufficient  information  to  permit  making  an 
evaluation . ) 

We  would  lixe  to  suggest  that  results  obtained  will  be  rather  in- 


sensitive to  the  choice  of  downstream  boundary  conditions.  There  are  two 
arguments : 

(1)  For  Tollmien-Schlichting  waves  the  group  velocity  (c  ) 
is  of  order  1/3U„.  Here  then  cg  ~ 500  cm/s.  Then  to 
travel  a distance  of  the  order  of  50  cm  from  the  ribbon 
to  the  downstream  edge  requires  ~ 1/10  s,  which  is 
approximately  15  periods  of  oscillations,  and  one  might 
not  even  follow  the  developments  much  further  in  time. 

(2)  The  essential  use  of  the  velocity  on  the  boundary  is  to 
compute  the  pressure  within  our  region.  This  can  be 
obtained  by  taking  a special  solution  of  the  Poisson 
equation  and  adding  to  it  a solution  of  Laplace's 
equation 


V2P  = 0 


chosen  to  satisfy  the  boundary  condition.  If  the  y 
variation  of  the  boundary  condition  is  ~ eiky,  then 
the  fall-off  in  the  x direction  is: 

-k| x-L| 

P ~ P e 1 
o 

Now  a reasonable  guess  for  the  y variation  is  k26  ~ 2 jt. 

(Remark:  We  are  thinking  of  deviations  from  the  Blasius 

profile.  For  the  latter,  of  course,  P is  independent 
of  y.)  This  suggests,  then,  that  the  boundary  conditions 
imposed  at  the  downstream  edge  will  have  little  effect 
a few  wavelengths  upstream.  It  is  for  this  reason  that 
we  have  suggested  that  the  downstream  edge  of  the  calcu- 
lation region  be  some  distance  beyond  where  we  hope  to 
find  a turbulent  transition. 

A further  argument  that  suggests  things  are  insensitive  to  the  down- 
stream boundary  conditions  is  found  by  noting  that  the  velocities  (phase 
and  group)  of  the  Tollmei n-Schlichting  waves  are  usually  small  (~  1/3) 
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compared  to  the  free-stream  velocity.  This  suggests  that  any  such 
disturbance  will  be  convected  downstream  and  have  little  effect  upstream. 

In  Appendix  C we  give  a simple  model  for  which  the  arguments  for 
exponential  fall-off  and  downstream  convection  are  verified  in  detail. 

However,  it  is  obvious  that  none  of  the  above  arguments  are  rigorous. 
It  is  essential  that  whatever  ooundary  conditions  are  used  should  be 
varied  (within  reasonable  limits)  to  see  how  far  upstream  we  need  to  be 
in  order  to  have  confidence  that  the  results  are  not  artifacts  caused  by 
the  assumptions  made. 

B . Calculations  at  Breakdown 

As  we  have  noted,  there  comes  a point  where  the  boundary  layer  is  so 
modified  by  the  nonlinear  wave  interactions  that  it  becomes  strongly  un- 
stable. At  this  point  high-frequency  oscillations  suddenly  appear.  In 
model  calculations  the  instabilities  grow  by  orders  of  magnitude  in  a 
fraction  (~  1/10)  of  the  primary  period.  The  scheme  described  above  will 
net  be  able  to  follow  such  rapid  variations.  At  best,  a large  irregular 
change  in  velocities  between  two  successive  time  steps  might  be  found. 
However,  even  this  may  not  occur — the  large  space  and  time  differences 
may  average  out  such  irregularities.  One  approach  to  proceeding  further 
is  suggested  by  the  fact  that  breakdown  (Klebanoff  et  al.,  1961)  appears 
to  occur  at  very  localized  points.  Hence,  we  can  istop  the  calculation 
when  at  some  x,z  there  is  a profile  that  according  to  linear  stability 
theory  is  highly  unstable.  Then  we  consider  only  a small  region  in  the 
vicinity  of  this  point.  For  this  region  we  introduce  a much  smaller 
space-time  mesh  and  then  integrate  the  Navier-Stokes  equation  with  initial 
conditions  obtained  from  the  previous  calculation. 
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Ill  COMPUTER  REQUIREMENTS 


A.  Introduction 

The  basic  problem  under  consideration  is  that  of  flow  over  a flat 
Plate  as  illustrated  in  Figure  6.  L,  M,  and  N are  the  numbers  of  mesh 


v (NORMAL  TO  PLATE), 
M-  32-64 


FIGURE  6 GEOMETRY  OF  FLAT-PLATE  NAVIER-STOKES  PROBLEM 

points  in  the  x,  y,  and  z directions  at  which  one  calculates  the  flow 
variables  v and  P.  Two  approaches  to  the  solution  of  the  Navier-Stokes 
equations  on  digital  computers  are  considered  here:  the  finite-difference 

scheme  of  Dr.  Chester  Grosch,  and  the  spectral  method  of  Dr.  Steve  Orszag. 
We  have  examined  the  finite-difference  method  in  considerable  detail 
with  the  help  of  Dr.  Grosch,  but  the  spectral  method  has  proved  more 
difficult  to  investigate,  as  noted  below.  The  objective  of  this  section 
is  to  make  some  computer-run-time  estimates  for  both  methods  on  three- 
dimensional  problems  of  interest.  The  results  should  not  be  taken  as  a 


comparison  of  the  two  methods  (as  we  discuss  below),  but  rather  as  an 
indication  of  the  computer  run  time  and  therefore  cost  that  one  can 
expect  on  optimized  "production  runs"  for  useful  problems. 

B.  Finite-Difference  Method 

The  finite-difference  method  for  flow  over  a flat  plate  in  two  di- 
mensions (x,y)  is  discussed  in  detail  by  Grosch  (1974).  The  method  is 

fairly  easily  generalized  to  three  dimensions.  The  basic  changes  are 

2 

that  the  Poisson  equation  for  the  pressure  V P-Q  = 0 must  be  solved  in 
three  dimensions  and  that  three  components  of  v are  required,  thus 
complicating  the  finite-difference  equations.  The  Poisson  equation  is 
solved  by  taking  a two-dimensional  Fourier  transform  in  the  (x,z)  plane 
for  each  of  the  M values  of  y in  the  mesh.  We  ther.  solve  (LN)  tri- 
diagonal equation  systems — i.e.,  one  system  for  each  mesh  point  in  the 
(x,z)  plane.  This  solution  yields  the  Fourier  components  of  P for  each 
plane  perpendicular  to  the  y axis.  An  inverse  Fourier  transfrom  then 
yields  P.  A set  of  three  finite-difference  equations  yields  the  three 
components  of  v and  the  loop  is  completed.  Table  1 summarizes  the 
operations  required  and  the  total  number  of  operations  per  time  step. 
Multiplying  the  terms  in  Table  1 by  their  appropriate  computation  times 

yields  the  total  computer  time  (AT  ) per  calculation  time  step: 

f d 

AT  * LMN  it  (95  + 3 log  LN)  + t (81  + 2 log  LN)  + 28t  1 (8) 

fd  ( a 2 m 2 s ) 

where  t , t , and  t are  the  particular  computer's  add,  multiply,  and 
am  s 

memory  access  times. 

C.  Spectral  Method 

It  is  impossible  to  determine  exactly  how  Orszag  implemented  his 
calculation.  If  his  report  (Orszag,  1974)  is  examined  it  will  be 
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THREE-DIMENSIONAL.  FINITE-DIFFERENCE  TECHNIQUE 


noticed  that  no  detailed  discussion  is  given  about  the  solution  technique. 
Verbal  conversations  with  Orszag  helped  some,  but  it  is  difficult  to 
judge  the  validity  and  efficiency  of  the  calculations.  In  at  least  one 
case,  the  report  contains  conflicting  statements  concerning  the  step- 
integration  method  employed.  "Adams-Bashforth"  is  indicated  on  p.  31, 
but  Runge-Kutta  is  used  in  the  actual  code.  Since  the  computer  code 
still  contains  bugs,'  it  is  not  clear  just  how  valuable  the  results  of 
the  study  were.  Although  it  is  not  evident  from  the  report,  it  is  abso- 
lutely true  that  Orszag  has  made  significant  contributions  to  the  field. 

The  general  method  can  be  described  as  a part  spectral  method  and 
part  finite-difference  method.  According  to  Orszag,  20  two-dimensional 
real  FFTs  (fast  Fourier  transforms)  and  two  one-dimensional  FFTs  are 
required  for  each  plane  normal  to  the  flow  per  time  step,  and  the  time 
spent  in  FFTs  is  60%  of  the  total  computer  run  time.  Also,  he  stated 
that  it  took  9 seconds  per  time  step  for  L = 128,  M = 8,  and  N = 64 . 

Thus,  a general  formula  for  estimating  the  run  time  on  a CDC-7600  can 
be  derived  as  follows. 

First,  let  us  consider  FFT  operations: 

(1 ) Complex  transform  in  one  dimension  requires 
j log^  (N)  butterflies  with 

4 multiplies  + 2 additions/rotation  + 4 additions/butterfly 
= multiplies  + 6 additions  for  each  butterfly  yielding 
3N  log^CN)  additions  and 
2N  log^(N)  multiplications  in  all. 

(2)  Real  transform  in  one  dimension  requires 
l0B2  (2)  + 1 butterflies  or 
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N[l°g2  (2)  + X]  “ N log2  N multiPlicationS  and 

(?) 


log2  N additions . 


(3\  2-p  complex  transform  (M  X N)  requires 


M 


£ los2  00 


followed  by 


M . 

— log  M 
2 2 


butterflies  (1-D  Transforms) 


butterflies  (1-D  Transforms) 


_ 

2 


log  M + log  N 
2 ^ 


= NM  log ^ MN  butterflies 


or 


(4  multiplication 
NM  log2  (MN)  X (e  addition 


(4  MN  log2  (MN)  multiplications 
= | 6 MN  log2  (MN)  addition. 

(4)  2-D  real  transform  (M  X N)  requires 


log  1 

(|) 4 ‘11 

I4 

[ 2 

\ 2/  J) 

(m 

/m\  J! 

111 

l0g2 

(?) * iii 

butterflies  followed  by 
butterflies 


= T [10g2  (2)  * l0g2  (1)  + 2]  ” ( 4 ) 


log2  MN  butterflies 


or 


, MN  log2  MN  multiplies 


iOK  mn  additions 
2 ?■ 


Orszag  has  estimated  that  the  (20L)  two-dimensional  real  FFTs  and 
(2L)  one-dimensional  FFTs  require  about  60%  ol  his  eo.puting  time.  For 
the  7600,  multiplication  Is  about  twice  the  time  of  an  addition.  Define 
one  operation  as  an  addition;  hence,  multiplication  Is  equivalent  to  two 
operations  and  a butterfly  Is  14  operations; 

23 


2L  (ID  FFTs)  * 14L  (“)  + (“)  log^ 


operations 


20L  (2D  FFTs) 


~ 20  X 14L 


log  (MN) 


operations 


The  total  run  time  (AT,)  per  computation  time  step  on  the  7600  is  then 


AT 


W R/  (o!e)  (14L)  5 m log2  (MN)  + log2M  f 1ob2N 

rs  * K'L  |ll7  UN  log2(MN)  + 5.8  M log  M + 5.8  N log  n| 


where  K is  a calibration  factor  (equal  to  the  add  time  in  seconds). 

Now,  from  Orszag's  experience  on  the  7600  we  know  that  AT  « 9S  for  L = 

s 

128,  M = 8,  and  N = 64.  Substituting  this  into  the  above  equation  yields 
K'  « 1.3  X 10  , so  the  effective  add  time  (t  ) is  about  one  ps . 

Now,  for  a single-calculation  time  step  we  will  allow  60%  of  the 
run  time  to  the  20L  two-dimensional  FFTs— i .e. , neglecting  the  two  one- 
dimensional FFTs,  we  allow  30%  of  the  run  time  to  the  solution  of  tri- 
diagonal equation  systems,  and  10%  of  the  run  time  to  14  memory  accesses 
per  mesh  point.  So  from  Orszag's  experience  we  can  derive  an  approxi- 
mate formula  giving  AT, , the  computing  time  required  for  a single  time 
step  in  the  Navier-Stokes  calculation: 


*TS  ~ LMNtg  (45  log2MN)  + LMNt  (30  log  MN)  + 14  LMNt  (9) 

M 6 S 

where  L,  M,  and  N are  the  numbers  of  mesh  points  in  the  x (flow), 

y (normal),  and  z (span)  directions;  and  t , t , and  t are  the  computer's 

a m s 

effective  addition,  multiplication,  and  memory  access  times. 
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D.  Comparison  of  Computing  Time  Estimates 


The  computing-time  estimates  (per  time  step)  for  the  finite  difference 

and  spectral  methods  discussed  above  are  given  by  Eqs.  (8)  and  (9).  The 

particular  computer's  effective  add,  multiply,  and  memory  access  times 

are  denoted  by  t , t , and  t . Two  features  are  common  to  both  methods: 
am  s 

(1)  The  computing  time  increases  approximately  linearly 
with  the  number  of  mesh  points. 

(2)  The  multiplication  terms  are  more  heavily  weighted 
than  the  addition  terms  when  one  considers  that 
multiply  times  are  typically  about  twice  as  long  as 
add  times. 

However,  there  are  also  some  interesting  differences.  The  finite- 
difference  method  is  more  sensitive  to  memory  access  time  and  therefore 
will  fare  relatively  worse  or.  a memory- limited  machine.  For  a typical 
application  the  spectral  method  is  most  heavily  governed  by  the  multi- 
plication time,  while  the  finite-difference  method  is  about  equally  de- 
pendent on  the  multiplication  and  memory  access  times. 

The  computing- time  estimates  given  below  should  be  qualified  on 

several  points.  First,  the  estimaces  are  intended  only  as  a very  rough 

guide  to  the  sort  of  computing  times  one  might  expect  on  rather  highly 

optimized  "production  runs."  We  have  used  effective  values  of  t^,  t^, 

and  t derived  from  the  experience  of  Dr.  Steve  Orszag  on  the  CDC  7600 
s 

computer  where  he  han  taken  care  to  reduce  the  run  time  by  optimization 
and  the  use  of  many  assembly-language-coded  subroutines.  The  finite- 
difference  method  of  . Chester  Grosch  has  never  been  optimized  for  the 
7600,  so  the  figures  below  are  in  a sense  a rough  projection  of  what  the 
finite-difference  method  could  do  if  optimized  in  a similar  fashion. 
Spectral  methods  are  thought  to  give  better  resolution  of  the  hydrodynamic 
variables  for  a given  mesh  size.  Thus,  for  calculations  of  a given 
accuracy,  the  number  of  mesh  points  required  could  be  less  for  the  spectral 
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method.  Finally,  the  algorithm  for  the  finite-difference  method  has  been 
examined  in  more  detail  than  has  that  for  the  spectral  method.  Given 
these  uncertainties,  one  should  regard  the  estimates  below  as  only  rough 
indications  of  the  range  of  computing  time  required,  and  not  as  a com- 
parison of  the  two  methods. 

For  problems  of  interest  the  numbers  of  mesh  points  required  in  the 

x,  y,  and  z directions  fall  in  the  following  ranges:  L ~ 128-256,  M ~ 

32-64,  and  N ~ 8-16.  In  Table  2 we  have  estimated  the  run  times  for 

maximum  (~  250,000  mesh  points),  minimum  (~  33,000),  and  typical  (~  66,000) 

problems.  AT,  the  computer  time  per  calculation  time  step  is  calculated 

from  Eqs . (8)  and  (9)  using  "effective"  values  of  t - 130  ns  t ~ 260  ns 

a m 

and  ts  ~ 1 ps.  These  effective  values  are  derived  from  the  experience 
of  Dr.  Steve  Orszag  on  the  CDC  7600  computer,  but  we  have  used  them  for 
estimating  both  AT_  and  ATfrf.  Typically,  300  to  500  time  steps  are  re- 
quired to  carry  a calculation  to  the  transition  stage,  so  we  have  taken 
T = 400  AT  as  the  total  run  time. 


Table  2 


ESTIMATED  COMPUTING  TIME  FOR  SPECTRAL 
AND  FINITE-DIFFERENCE  METHODS 


Problem 

L 

M 

N 

LMN 

AT 

fd 

AT 

s 

T 

fd 

T 

s 

Maximum 

256 

64 

16 

262,144 

20  s 

40  s 

2.0  hr 

4.0  hr 

Minimum 

128 

32 

8 

32,768 

2 s 

4 s 

0.2  hr 

0.5  hr 

Typical 

128 

64 

8 

65,536 

5 s 

8 s 

0.5  hr 

0.9  hr 

Note:  L - Along  flow;  M = Normal  to  flow;  N = Along  span 


26 


Again  .0  should  emphasize  that  those  figures  should  not  he  taken  as 
a comparison  of  the  finite-difference  and  spectral  methods- they  are  too 
rough  The  main  point  here  is  that  for  suitably  optimized  codes  the  run 
times  are  not  unacceptably  long.  Given  a typical  cost  of  - $700  per  hour 
for  a 7600,  the  production  run  cost  for  the  -typical''  problem  is  about 
$500.  Optimization  and  assembly  language  coding  (assumed  here),  which 
take  advantage  of  the  architecture  of  a given  computer,  can  give  a time 
reduction  as  large  as  a factor  of  five.  So  one  could  ezpect  in  unopti- 
mlzed  run  to  cos,  about  two  to  five  times  as  much  as  these  "production 


runs 
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IV  WAKE  OF  A FLAT  PLATE 


A prob  em  closely  related  to  the  above  transition  in  the  boundary 
layer  of  a flat  plate  is  that  of  transition  in  the  wake  of  a flat  plate. 
Some  reasons  for  considering  this  problem  next  are: 

(1)  The  calculations  are  very  similar  to  the  previous  one. 

A suitable  program  for  the  flat-plate  calculation 
would  require  modest  changes  in  order  to  treat  the  wake. 

(2)  Detailed  experimental  information  is  available  for 
comparison  (for  example,  see  Sato  and  Kuriki  , 1961). 

(3)  The  nature  of  the  transition  is  rather  gentle.  In 
contrast  to  many  other  turbulent  transitions,  this 
one  is  not  achieved  by  the  development  of  sharp 
bursts.  Accordingly,  it  may  be  expected  that 
transition  can  be  followed  numericr  lly  very  far. 

Drawbacks  to  this  calculation  are: 

(1)  As  previously,  the  downstream  boundary  conditions 
are  ambiguous . 

(2)  The  applications  to  practical  design  are  somewhat 
indirect . 

Our  understanding  of  this  transition  runs  as  follows:  The  Goldstein 

(laminar)  wake  (Goldstein,  1933)  is  unstable  according  to  two-dimensional 
linear  stability  theory.  Thus  a small  perturbation  is  amplified  when 
one  proceeds  downstream.  Eventually  nonlinear  effects  take  over.  (The 
growth  rate  is  not  the  predicted  exponential,  and  harmonics  of  the  in- 
duced perturbations  appear.)  Further  downstream  a distinct  three- 
dimensional  pattern  appears  and  the  flow  becomes  more  and  more  irregular. 
However,  we  emphasize  that  no  sharp  bursts  or  spikes  appear.  This  should 
be  an  ideal  problem  for  numerical  simulation. 


Preceding  page  blank 
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Since,  as  indicated,  the  problem  is  so  much  like  the  flat-plate 
transition , computing  requirements  will  also  be  pretty  much  the  same 


V WEDGE  FLOWS 


I 

J 

The  numerical  simulation  of  transition  in  wedge  flows  would  be  of 
considerable  interest.  The  reasons  for  this  are: 

(1)  The  calculation  is  very  similar  to  the  flat  plate. 

■.Indeed,  one  code  (Grosch,  1974)  has  been  written 
to  include  this  possibility.]  The  essential  reason 
for  this  is  that  the  boundary-layer  equations  for 
the  flat  plate  are  just  a special  case  of  that  for 
the  wedge.  Thus,  if  pjt  is  the  included  angle,  the 
potential  flow  is 

U(x)  = UQxm  , UO) 

where 


m = P/(2  - B) 


and  the  boundary  layer  is  described  by 
(Schlichting,  I960,  p.  143): 

the  equations 

m § 

U x f (Tl) 
o 

|m  + 1 ..  m-1 

V‘-<l  2 VV 

If  + m-  1 f'| 

I m + 1 ) 

» 

(ID 

T]  = 

. u 

m + 1 o m-1 

(12) 

2 v ^ 2 

f + ff"  + p[i  - (f  ')2]  = o 

1 

(13) 

subject  to  f(0)  = f'(0)  = 0 , f'(co)  = 1 

We  note  that  with  0 = 0 this  is  just  our  flat-plate 
equation,  (The  factor-of-2  difference  occurring 
between  Eqs . (5)  and  (13)  is  due  to  a different 
normalization.) 


{ 


1 


(2)  By  varying  the  parameter  (3  we  can  study  the  effects  of 
pressure  gradient.  Thus,  with  p > 0 we  have  accelerated 
flow  (and  therefore  greater  stability),  while  for  P < 0 
we  have  deceleration  (and  more  instability). 

(3)  With  P = 1/2  we  have  the  equations  of  the  boundary  layer 
for  a rotat tonally  symmetric  flow. 

(4)  Measurements  of  the  boundary- layer  transition  on  a flat 
plate  with  mild,  favorable  pressure  gradients  (DeMetz, 
1974)  are  available  for  comparison  with  the  calculations. 

As  we  have  indicated,  the  computing  requirements  are  as  for  the 


flat-plate  case . 


VI  POSSIBLE  FUTURE  PROBLEMS 


There  seem  to  be  many  different  mechanisms  for  transition  to  turbu- 
lence. For  an  understanding  of  these  mechanisms,  it  would  be  desirable 
to  do  numerical  simulation  in  relatively  pure  situations.  Here  we  give 
a list  of  some  of  these  and  an  indication  of  why  they  are  of  interest. 

(Since  the  calculations  are  of  a somewhat  different  nature  than  those 
described  above  and  have  less  direct  practical  application,  we  think 
they  should  be  done  at  some  later  date.  Accordingly,  we  have  not  made 
any  computing  estimates.  However,  it  seems  that  they  are  problems  of  a 
magnitude  similar  to  that  of  the  problems  discussed  above~but  computa- 
tion would  probably  require  new  codes.) 

A.  Couette  Flow 

This  flow  between  concentric  rotating  cylinders  is,  of  course, 
classic.  Extremely  good  experimental  results  are  available  (Donnelly, 

1963;  Coles,  1965).  From  a theoretical  standpoint  this  is  a very  clean 
problem.  The  ambiguity  of  the  boundary  conditions  disappears.  In  the 
radial  direction  we  have  well  defined  boundary  conditions  on  the  cylinders. 
Azimuthally  we  have  periodicity.  In  the  longitudinal  direction  periodicity 
is  a reasonable  requirement. 

When  the  inner  cylinder  rotates  more  rapidly,  there  is  seen  experi- 
mentally a rich  set  of  phenomena  (Coles,  1965).  When  the  rotation  rates 
are  changed  we  pass  from  the  original  laminar  state  through  a multitude 
of  different  states  with  well  organized  flow  patterns.  In  addition, 
hysteresis  effects  manifest  themselves--i .e . , the  state  achieved  depends 
on  the  past  history.  It  would  be  fascinating  to  reproduce  these  effects 
by  numerical  simulation  and  it  should  be  possible  to  do  so. 


When  the  outer  cylinder  rotates  more  rapidly,  something  that  Coles 
(1965)  has  called  "catastrophic  transition"  occurs.  There  is  a sudden 
jump  to  states  with  significant  turbulence.  Frequently  there  are  alternate 
spirals  of  turbulent  and  laminar  fluid.  This  could  well  be  not  resolvable 
with  present  computers. 

B.  The  Benard  Problem 

This  is  the  problem  of  a fluid  heated  from  below.  Mathematically 
it  is  closely  related  to  the  Couette  f.ow.  Linear  stability  theory 
correctly  predicts  the  change  from  the  quiescent  state.  However,  it 
makes  no  prediction  as  to  the  structure  that  develops — hexagonal  cells, 
rolls,  etc.  Model-type  calculations — essentially  including  some  non- 
linearity— have  been  done  to  determine  the  structure.  Since  these  cal- 
culations are  quite  similar  to  those  made  to  describe  the  nonlinear 
region  in  transitions  on  a flat  plate,  it  would  be  worthwhile  to  see 
with  numerical  simulations  how  good  these  calculations  are. 

C.  Poiseuille  Flow 

1 . Between  Parallel  Planes 

The  results  of  linear  stability  theory  (Figure  7)  for  this  flow 
are  strikingly  similar  to  those  for  tne  flow  over  the  flat  plate.  Accord- 
ingly, it  may  be  expected  that  transition  to  turbulence  will  occur  as  in 
the  case  of  the  flat  plate.  One  major  advantage  is  that  the  downstream 
boundary- condition  problem  can  be  avoided  by  the  reasonable  assumption 
of  periodicity.  A disadvantage  compared  to  the  flat-plate  case  is  that 
comparable  detailed  experimental  information  does  not  seem  to  be  available. 
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FIGURE  7 THE  FORM  OF  THE  NEUTRAL  CURVE  IN  THE  (k,Re) 
PLANE  FOR  A PLANE  POISEUILLE  FLOW.  Source: 
Monin  and  Yaglom  (1971). 


2.  In  Pipes 


This  is  a particularly  interesting  situation  since,  in  contrast 
to  all  other  problems  we  have  discussed,  the  flow  is  apparently  stable 
with  respect  to  infinitesimal  disturbances.  We  say  apparently  since: 

(1)  Extensive  theoretical  efforts  using  linear- 
stability  theory  have  failed  to  find  any 
unstable  modes. 

(2)  Experimentally  the  transition  Reynolds'  number 
seems  to  go  ever  higher,  the  lower  the  ambient 
disturbance  level.  (Of  course,  neither  of  these 
is  a proof.) 

At  least  two  ideas  as  to  how  the  observed  turbulence  originates 
have  been  suggested.  One  is  that  there  is  a finite  amplitude  instability. 
If  this  is  the  case,  numerical  simulation  seems  to  be  of  little  value. 

The  observed  transition  in  terms  of  turbulent  bursts  suggests  that  in  a 
calculation  we  would  either  be  in  a laminar  region  (of  little  interest), 
or  in  a turbulent  region  where  we  cannot  resolve  the  flow.  The  second 
suggestion  (Smith.  1960)  is  that  turbulence  originates  in  the  inlet 
region  before  the  full  Poiseuille  profile  develops.  If  this  picture  is 
correct,  the  numerical  calculations  would  be  similar  to  those  for  the 
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flat-plate  boundary  layer  It  stmi.u  k 

Sh0l'ld  be  Possible  to  calculate  most  of 

the  .ay  to  transition  to  fully  developed  turbulence. 

“•  Modifications  of  the  Flat-Plate  Prohl— 

Assuming  successful  completion  of  the  flat-plate  calculations.  „ 
number  interesting  numerical  experiments  might  be  performed.  Some 
require  very  little  code  modification.  Thus,  the  effects  of  suction 
toughness,  and  compliant  boundaries  can  be  treated  by  merely  changing 
c oundary  conditions  imposed  on  the  plate.  The  effects  due  to  heating 
could  also  be  treated.  An  additional  equation  <for  the  temperature, 

•ould  have  to  be  used,  and  the  appropriate  variation  in  viscosity  needs 
O deluded.  While  some  recoding  would  be  necessary,  the  order  of 
magnitude  of  the  computation  needed  would  not  he  changed. 
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VII  ON  CONTROLLING  THE  TRANSITION  TO  TURBULENCE 


* 


( 


The  development  of  a three-dimensional  numerical  simulation  code 
will  give  us  a powerful  new  tool  for  assessing  methods  to  control  the 
transition  to  turbulence  on  a wide  range  of  bodies— i.e.,  not  Just  bodies 
of  revolution  and  two-dimensional  airfoils.  We  use  the  phrase  "controlling 
the  transition"  in  a general  sense  so  that  studies  of  effects  that  sub- 
stantially advance  the  transition  (e.g. , excessive  surface  roughness) 
are  also  included  in  the  scope  of  the  investigations. 

Let  us  first  turn  to  ways  to  calculate  and  control  linear  Tollmein- 
Schlichting  instabilities.  Here  a number  of  methods  have  already  been 
suggested:  (1)  compliant  surfaces,  (2)  changing  the  equation  at  state 

(e.g.,  adding  heat,  polymers,  etc.),  (3)  suction,  and  (4)  control  of  the 
pressure  gradient.  Except  for  pressure  gradients,  which  have  enjoyed 
wide  use,  the  optimum  ways  to  use  Methods  1 through  3 remain  to  be  devised. 
Furthermore,  even  though  one  is  dealing  with  linear-stability  problems, 
there  are  important  nonlinear  effects  in  the  ways  the  laminar  boundary 
layer  and  potential  flow  patterns  are  altered  by  the  measures  used  to 
delay  the  transition.  For  example,  adding  heat  changes  the  viscosity, 
and  the  viscosity  gradient  in  turn  alters  the  laminar  boundary  layer  in 
a direction  toward  stability.  Clearly,  heat  added  at  one  location  on  the 
surface  of  a body  will  diffuse  across  the  boundary  layer  and  have  a much 
smaller  effect  downstream.  Thus  the  combination  of  heat  transport  via 
diffusion  and  convection,  coupled  with  the  linear-instability  mechanisms, 
becomes  a complicated  problem  that  is  probably  most  easily  solved  by  a 
direct  simulation  rather  than  linear-stability  analysis.  The  present  codes 
must  be  generalized  to  include  heat  transport,  of  course. 


I 


3 


4 


A 
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The  second-class  methods  to  alter  the  transition  involve  nonlinear 
changes  in  the  flow  within  the  boundary  layer — e.g.,  nonparallel  flows, 
surface  roughness,  and  abrupt  changes  in  the  body  shape.  As  examples 
we  can  cite  the  early  transition  to  turbulence  that  occurs  on  some  rough, 
blunt  bodies  and  the  use  of  small  span-wise  or  flow-wise  grooves  to 
stabilize  the  Tollmein-Schlichti ng  waves.  Also,  swept-wing  airfoils 
have  boundary  layers  with  nonparallel  flow.  The  numerical  simulation 
approach  gives  us  a single  technique  that  can  handle  this  wide  variety 
of  problems. 

Overall,  our  studies  indicate  that  numerical  simulation  techniques 
have  sufficient  resolution  to  accurately  predict  transition  and  that  they 
can  be  adopted  to  circumstances  where  the  laminar  boundary  layers  differ 
appreciably  from  the  simple  parallel-flow  ones  now  commonly  used. 

Initial  studies  of  such  phenomena  would  take  a year  to  complete,  and  a 
comprehensive  program  would  require  three  to  five  years. 
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VIII  THE  IMPACT  OF  NUMERICAL  SIMULATION  OF  TURBULENCE 


Let  us  suppose  that  three-dimensional  numerical  simulation  technique 
have  successfully  reproduced  many  of  the  phenomena  that  occur  in  flat- 
plate  flows.  What  comprises  the  next  generation  of  problems  both  in  the 
practical  and  research  areas? 

An  important  class  of  practical  problems  concerns  the  simulation  of 
the  transition  to  turbulence  in  configurations  where  linear-stability 
analysis  fails  or  is  excessively  complicated  due  to  geometry.  These 
problems  will  require  the  development  of  a more  flexible  code  that  can 
handle  both  fairly  general  boundary  conditions  (e.g.  , surfaces  rough  on 
the  scale  of  the  boundary  layer)  and  heat  transport.  Fortunately,  almost 
every  body  can  be  represented  by  a sequence  of  simulations  each  of  which 
represents  a small  portion  of  the  surface. 

The  impact  of  the  numerical  simulation  of  turbulence  is  that  it  will 
be  a single  technique  that  will  handle  many  problems  not  only  of  a linear 
stability  type  but  also  inherently  nonlinear,  such  as  finite-amplitude 
surface  roughness.  The  ability  to  compute  the  combined  effects  of  both 
linear  and  nonlinear  processes  for  initiating  transition  is  a strength 
unique  to  the  numerical  simulation  approach. 

The  rext  generation  of  research  problems  should  be  chosen  to  help 
us  understand  how  well  the  simulation  code  represents  fully  developed 
turbulence.  Estimates  clearly  indicate  that  the  inertial  range  will  be 
described  poorly  if  at  all.  For  example,  high-frequency  bursts  of 
turbulence  may  not  occur  in  the  correct  frequency  range.  On  the  other 
hand,  many  of  the  overall  properties  of  a turbulent  boundary  layer  (such 
as  drag,  heat  transport,  etc.)  may  be  quite  well  represented  since  they 


depend  much  more  on  the  larger,  energy-containing  eddies.  Therefore, 
there  exists  a strong  possibility  that  numerical  simulation  may  produce 
a useful  enough  representation  of  the  turbulent  boundary  layer  so  that 
the  effect  of  body  shape,  etc.,  on  the  transport  properties  of  the  layer 
can  be  calculated — at  least  in  a semi-quantitative  manner.  Such  a pro- 
gram would  require  codes  with  flexible  boundary  conditions  and  carefully 
thought  out  sub-grid  closure  schemes.  Its  impact  would  be  to  produce  a 
first-principles  calculation  of  the  arbitrary  coefficients  that  occur  in 
turbulent-boundary-li yer  theory,  and,  more  importantly,  to  see  if  there 
are  circumstances  in  which  these  coefficients  vary  by  a significant 
amount . 


IX  RECOMMENDATIONS 


A program  of  numerical  simulation  that  could  be  expected  to  have  a 
significant  impact  on  our  knowledge  of,  and  ability  to  predict,  the 
transition  from  laminar  to  turbulent  flow  might  consist  of  the  following: 

(1)  A commitment  to  the  program  for  a significant  length  of 
time.  Crash  programs  seem  to  be  counterproductive.  Four 
to  five  years  is  a reasonable  time  scale. 

(2)  An  aspect  of  competition.  There  should  be  at  least 
two  groups  performing  the  calculations.  In  the  beginning 
they  should  be  calculating  as  closely  as  possible  the 
same  problem,  and  with  the  same  parameters  used  in  the 
experiments.  Later,  when  confidence  has  been  obtained 
the  groups  could  begin  dividing  up  the  problems  suggested. 

It  is  reasonable  that  the  groups  use  somewhat  different 
computation  schemes.  A number  of  such  schemes  suggest 
themselves.  It  would  be  efficient  to  use  the  program 
to  evaluate  their  relative  merits.  The  schemes  used  to 
date  i.e.,  those  of  Grosch  and  Orszag — seem  sufficiently 
representative. 

(3)  A small  advisory  panel.  This  might  consist  of  a group 
(~  4)  of  active  workers  in  the  field.  The  members  might 
be  an  experimentalist,  a theoretical  hydrodynamicist , a 

l#  computer  expert,  and  a numerical  analyst.  They  would 

be  expected  to  work  closely  with  the  computing  groups, 
Hopefully,  they  would  spot  difficulties  or  uncertainties 
* in  given  calculations  and  direct  the  program  (i.e.,  while 

we  have  given  a tentative  ordering  of  problems  to  be 
done,  it  would  be  expected  that  the  advisors  would  con- 
| tinually  reassess  the  next  stages  of  the  program  in  the 

light  of  what  had  been  accomplished). 

(4)  An  expenditure  of  the  order  of  $200,000  per  year.  This 
might  produce  between  30  and  100  runs  on  one  of  the 
problems  outlined  above.  We  envisage  a few  man-years 

effort  and  computation  on  a computer  comparable  to 
the  7600.  This  sum  may  or  may  not  include  support  for 
experimental  work  that  might  be  found  advisable. 


41 


Appendix  A 

LINEAR-STABILITY  ANALYSIS  AS  A PRACTICAL  TOOL 


Preceding  page  blank 


43 


Appendix  A 


LI  NEAR- STABILITY  ANALYSIS  AS  A PRACTICAL  TOOL 

On  a properly  designed  and  fabricated  body,  the  transition  from 
laminar  to  turbulent  boundary  layers  occurs  through  linear  instabilities 
whereby  free-stream  fluctuations  in  velocity,  pressure,  and  temperature 
enter  the  boundary  layer  and  are  spatially  amplified  via  the  Tollmein- 
Schlichting  waves.  A sufficiently  large  amplification  results  in  non- 
linear effects  that  initiate  the  transition  to  turbulence.  The  transition 
prediction  method  of  A.M.O.  Smith  and  his  colleagues  is  the  only  predic- 
tion method  that  now  makes  use  of  the  spatial  amplification  concept— a 
process  that  must  occur  physically.  Indeed,  Smith's  e method  has  been 
the  most  successful  predictor  in  ARPA's  program  as  well  as  in  predicting 
transition  for  the  various  airfoils. 

But  Smith's  method  has  not  received  wide  use  because  the  spatial 
amplification  rates  are  derived  from  the  solution  to  an  eigenvalue  problem 
involving  4th-order  differential  equations.  The  solution  to  these  equa- 
tions, which  must  be  done  by  computer  for  an  arbitrary  laminar  boundary 
layer,  are  conceptually  straightforward  but  quite  laborious  in  practice. 

A well  documented  and  widely  available  computer  program  is  called  for. 

We  therefore  recommend  that  ARPA  fund  an  effort  to  produce  a computer 
program  (most  likely  in  FORTRAN)  that  will  run  on  a number  of  the  more 
generally  available  scientific  computers.  The  core  of  the  program  will 
be  the  4th-order  eigenvalue  solver.  Auxiliary  parts  will  integrate  the 
spatial  amplification  factor  along  the  body,  and  calculate  the  laminar 
boundary  layer,  potential  flow,  etc. 
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The  wide  availability  of  such  a program  would  allow  body  designers 
to  test  proposed  bodies  against  Smith's  criterion  on  an  in-house  basis, 
hopefully  allowing  a more  rapid  convergence  to  optimum  designs. 

The  need  for  this  computer  program  is  an  example  of  a situation  in 
which  an  increased  understanding  of  the  basic  physics  has  led  to  more 
sophisticated  design  criteria  that  can  be  implemented  only  by  a computer 
program.  Of  course,  some  generality  is  lost  from  the  older,  more  gross 
engineering  criteria.  But,  it  is  only  through  such  a program  that  the 
effects  of  the  details  of  body  shape  on  the  transition  to  turbulence  can 
be  quickly  evaluated. 

In  closing  this  section,  we  should  remind  the  reader  that  transition 
can  occur  through  channels  other  than  linear  Tollmein-Scnlichtir.g  ampli- 
fication. High  levels  of  free-stream  turbulence,  excessive  surface 
roughness , and  poor  body  fabrication  are  examples.  The  spatial  amplifi- 
cation method  often  does  not  work  in  such  cases.  But  it  does  have  a good 
record  of  success  in  those  projects  where  designers  have  taken  the 
trouble  to  eliminate  these  other  sources  of  turbulence. 
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REMARKS  ON  A SPECIALIZED  NAVIER- STOKES  COMPUTER 
1.  Introduction 

Reference  2 has  considered  the  question  of  what  can  be  accomplished 
with  modern  machines  such  as  the  ILLIAC  IV,  IBM-370/195,  and  CDC-7600. 

The  purpose  of  this  appendix  is  to  consider  the  wisdom  of  designing  a 
special  computer  that  could  greatly  exceed  the  capabilities  of  the  current 
machines  in  solving  only  a specialized  class  of  problems — i e,,  the 
Navier-Stokes  equations,  and  similar  partial  differential  equations 
arising  in  plasma  physics,  hydrodynamics,  weather  prediction,  geology, 
etc. 


a . An  Interesting  Calculation 

An  interesting  numerical  simulation  well  beyond  the  capabilities 

of  current  machines  is  the  fully  resolved  three-dimensional  turbulent 

4 

fluid  flow  for  a Reynolds  number  of  10  . This  requires  a grid  of 
4 9/4  9 

(10  ) =10  mesh  points  for  each  point  in  time.  Each  mesh  point  must 

represent  the  velocity  and  pressure  fields  so  that  three  velocities  and 
one  pressure  term  are  required  at  least  for  the  present  instant  in  time 
and  the  immediate  past.  The  past  values  are  needed  in  order  to  march 

Q 

forward  in  time  with  a numerical  integration  procedure.  Thus  8 X 10 

scalar  variables  are  required,  and  about  (104)3  In  104  ==  1013  arithmetic 

operations  are  required.  A more  detailed  count  of  arithmetic  operations 

of  an  efficient  finite  differencing  scheme  for  10  grid  points  indicated 

that  5 X io  additions  and  3 X 10  multiplications  are  needed  per  time 

12 

step.  A rough  guess  ic  that  about  10  instructions  would  be  required 
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n 


12 

to  implement  these  operations  so  that  about  2 X 10  read  accesses  from 
storage  are  required  for  each  simulated  time  step.  These  parameters  are 
summarized  in  Table  B-l. 


Table  B-l 

PARAMETERS  OF  AN  INTERESTING  TURBULENCE  CALCULATION 


Three-Dimensional  Geometry 
Reynolds  number: 

4 

R = 10 

Resolution: 

10®  spatial  mesh  points  per  time 

step 

Memory  size: 

1010  scaler  variables 

Arithmetic  operations 
(finite-difference  method) 

Additions: 

5 X 1011  per  time  step 

Multiplications: 

3 X 1011  per  time  step 

Instruction  execution  rate: 

1012  per  time  step 

Memory  bandwidth: 

2 X 1012  read  accesses  per  time 

step 

b , Current  General  Purpose  Computers 

Table  B-2  lists  some  critical  state-of-the-art  parameters  of 
modern  computers.  Note  that  with  a simple  array  system  of  secondary- 
storage  access  systems,  more  than  2 X 10  X (2  X io  ) =10  s would 

be  required  Just  to  update  each  grid  point  for  each  simulated  time  step. 
Parallel  access  to  the  secondary  storage  will  reduce  this  time,  however. 

The  instruction  execution  time  for  a very-high-performance  single  CPU 

-9  9 w 

(control  processing  unit)  is  at  best  25  X 10  s,  so  that  25  X 10  X 

2 X 1012  = 5 X io4  s or  about  14  hours  per  simulated  time  step  would  be 

9 

required.  Thus  it  seems  infeasible  to  simulate  10  grid  points  since 
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Table  B-2 


MAXIMUM-PERFORMANCE  PARAMETERS 
FOR  A IODERN  SINGLE,  PIPELINED  AND  CACHE  CPU  COMPUTER 


Direct-access  secondary-storage  capacity: 
Direct-access  secondary-storage  b:  idwidth: 
Arithmetic  operation  time: 

Cache  memory  access  time: 


1011  bits  maximum — i.e., 
about  2 X 10®  variables 

107  bits/s  maximum — i.e., 
about  2 X 106  variables/s 

25  X 10“9  s minimum  (pipe- 
lined) 

25  X 10-9  s minimum 


many  time  steps  would  be  required  in  a practical  problem.  Clearly,  the 
usual  von  Neumann  computer  architecture  cannot  be  used,  and  other  choices 
should  be  considered. 

c.  ILLIAC  IV 

The  ILLIAC  IV  (Bell  and  Newell,  1671)  is  a radical  departure 
from  classical  computer  archil  itn re.  It  was  designed  for  array- type 
calculations.  The  critical  parameters  of  the  ILLIAC  IV  are  given  in 
ble  B-3.  The  most  limiting  factor  is  the  limited  secondary  storage  of 
only  10^  bits  as  compared  to  the  10  bits  that  are  required.  Even  if 
the  ILLIAC  should  be  suprlemented  with  large  disc  drives,  the  bandwidth 
to  secondary  storage  would  be  less  than  10°  bits/s  as  compared  to  60  X 10 
bits  needed  for  one  simulated  time  step,.  It  would  take  more  thar  60  X 
10'9  = 600  s just  to  read  in  data  for  one  simulated  time  step.  Calcula- 
tions would  require  about  2 X 10  X 2 X 240  X 10  + 64  = 15,000  s per 

simulated  time  step,  if  all  input-output  and  grid-point  data  interchanges 
are  ignored.  Thus  even  an  extensively  enhanced  ILLIAC  IV  cannot  ue  us 
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Table  B-3 


ILLIAC  IV  PARAMETERS 


64  processing  elements 

Add  time  (each  element) : 

Multiply  time  (each  element): 

2048-64  words  of  240-ns 
memory /e lemen  t 

Secondary  storage  capacity: 

Secondary  storage  bandwidth: 


240  ns 
400  ns 


109  bits 
q 

10  bits/s  (to  all  elements  In  parallel) 


to  attack  this  problem,  and  other  types  of  computer  structures  need  to 
be  examined. 

2.  Computer  Architecture  for  Array  Calculations 
a.  Solution  Methods 

Nonlinear  partial  differential  equations  can  be  solved  by  a 
variety  of  techniques  such  as  finite  difference,  spectral,  psuedospectre 1 , 
and  Monte  Carlo  methods.  A special  computer  designed  for  any  one  of  these 
methods  will  not  be  optimal  for  any  of  the  others.  Also,  while  it  appeals 
that  finite-difference  methods  are  the  most  flexible  relative  to  boundary 
complexity,  etc.,  the  pseudospectral  methods  are  much  better  when  regular 
surface  boundaries  exist.  Furthermore,  integration  algorithms  nre  in  a 
dynamic  state  of  development  and  better  algorithms  are  quite  likely  to 
emerge  in  the  near  future.  Thus,  only  very  flexible  and  general  computer 
structures  should  be  considered. 
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b.  Computing  Requirements 


What  are  the  common  characteristics  of  a large  class  of  non- 
linear partial  differential  equation  solution  methods?  First  of  all, 
for  very  arbitrary  boundary  conditions  some  finite-difference  and  suc- 
cessive over-re xaxation  methods  must  be  employed.  However,  if  well 
behaved  boundaries  (planes,  spheres,  cylinders,  etc.)  exist,  then  the 

pseudospectral  methods  offer  large  savings  in  computing  effort.  In  fact 

g 

it  appears  that  the  10  grid  problem  will  be  attacked  in  the  foreseeable 
future  only  by  pseudospectral  or  other  methods  that  „ake  advantage  of 
particular  features  of  a given  problem.  Finite-difference  and  relaxation 
methods  require  only  that  information  be  passed  between  neighboring  grid 
points.  The  spectral  techniques  on  the  other  hand  need  global  informa- 
tion at  every  grid  point  in  order  to  calculate  a Fourier  or  other  trans- 
form (Orszag,  1974),  However,  a very  large  class  of  transforms  can  be 
implemented  through  the  fast  Fourier  transform  (FFT) . Thus  the  archi- 
tecture of  the  computer  should  be  such  that  grid  points  can  interchange 
information  with  their  close  neighbors  and  also  be  able  to  perform  3- 
dimensional  FFTs  on  the  quantities  of  the  grid  points. 

F~om  Table  B-l  it  is  evident  that  to  solve  the  Navier-Stokes 
equations,  for  a large  number  of  grid  points,  about  ten  variables  per 
grid  and  1000  operations  per  grid  point  are  required  for  each  simulated 
time  step.  More  complicated  problems  involving  compressible  fluids, 
electromagnetic  fields,  etc.,  might  double  both  of  the  quantities  so  that 
in  designing  a relati  'ely  flexible  computer,  perhaps  32  variables  per  grid 
point  should  be  selected  and  2000  to  4000  operations  per  grid  point  per 
time  step  assumed.  If  it  is  desired  that  a time  step  be  simulated  in 
about  a second,  then  about  3200/32  = 100  operations  per  second  are  needed 
for  every  variable. 


c. 


Performance  Tradeoffs 


Very  fast  pipelined  logic  (a  high-performance  processor)  could 
perform  about  2 X 10?  operations  per  second  while  a slow  serial  processor 
might  only  be  able  to  do  about  3 X lO3  operations  per  second. 

For  one  second  of  execution  time  per  simulated  time  step, 

2 X 10  /2000  — 10  grid  points  could  be  processed  by  a single  high- 

4 

performance  processor  with  (10  to  32)  X 10  registers  per  processor 

0 

(i.e.,  approximate  10  bytes).  However,  an  interconnected  array  of 

x 10^  or  10^  X 10^  high-performance  processors  would  be  required  to 
handle  the  109  grid-point  problems.  This  approach  seems  to  be  far  beyond 
near-term  technology  hence  is  therefore  considered  infeasible.  Even  if 
the  calculation  time  were  increased  considerably  and  the  logic  speeds 
could  be  greatly  improved,  an  impossibly  large  array  of  very-high- 
performance  computers  would  still  be  required,  in  any  case,  the  pro- 
cessors of  such  an  array  would  be  very  similar  in  structure  tc  the  present 
high-performance  machines  discussed  earlier  so  that  beyond  creating 
special  FFT  hardware  attachments  to  the  usual  sort  of  pipelined  processor 
(or  array  processors  for  ILLIAC  IV) , it  appears  that  little  can  be  gained 
by  a special  computer  architecture  using  high-performance  operations. 

An  interesting  situation  develops,  if  however,  the  low- 
performance  serial  organization  is  considered.  If  one  slow  serial  pro- 
cessor were  "devoted"  to  one  grid  cell  and  hence  32  variables,  the  re- 
sulting computer  structure  would  be  a three-dimensional  array  of  slow 
serial  processors,  each  with  about  32  registers,  connections  to  its 
neighbors  (six  each),  and  connections  to  implement  an  FFT  (see  Appendix 
E for  duration  of  the  FFT)  in  three  dimensions  (six  additional  connections 
per  cell). 

A more  efficient  arrangement  in  terms  of  total  number  of  inter- 

3 

connections  is  to  use  one  processor  for  every  (N  ) grid  points;  for  the 
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lowest-speed  serial  processor,  this  could  require  as  long  as  ten  seconds 
to  simulate  one  time  step.  This  may  very  well  be  an  optimal  arrangement 
for  a special-purpose  computer  devoted  to  solving  hydrodynamic  equations. 
As  N is  increased,  the  arrangement  of  cells  is  not  changed  but  a few 
more  instructions  per  cel1,  are  required  for  individual  grid-point  selec- 
tion within  a given  cell.  Table  B-4  illustrates  the  possibilities  for 
several  values  of  N.  It  is  significant  to  note  that  by  1980,  it  may  be 
feasible  to  construct  a computer  to  solve  the  10  grid  problem  for  about 
the  cost  of  the  ILLIAC  IV  (in  1970  dollars). 


3.  Conclusions 

It  can  be  concluded  that: 

(1)  Direc*-  numerical  simulation  of  turbulence  for  most 
important  and  interesting  problems  cannot  be  accomplished 
with  present-day  computers  or  even  with  general-purpose 
computers  proposed  for  the  near  future.  Hence,  many  of 
these  simulations  will  be  impossible  without  soecial- 
purpose  computers  that  have  been  custom  designed  to 
solve  partial  differential  equations. 

(2)  Algorithms  for  solving  partial  differential  equations 
are  in  a dynamic  state  of  development.  Thus  any  pro- 
posed special-purpose  computer  must  be  very  flexible 
relative  to  solution  method  as  well  as  to  the  kinds  of 
problems  it  can  handle.  In  particular,  it  should  be 
capable  of  finite-dif fere  ice  calculations  and  relaxation 
calculations  as  well  as  spectral  and  pseudo  spectral 
calculations  for  sets  of  partial  differential  equations 
that  arise  in  hydrodynamics,  plasma  physics,  aerodynamics, 
global  weather,  etc.  Hence,  each  grid  point  must  be 

able  to  exchange  information  with  neighboring  grid 
points  and  must  be  capable  of  in-place,  three-dimensional 
fast  Fourier  transformation  of  the  grid-point  quantities. 

(3)  A special  purpose  very-high-performance  computer  could 
be  organized  for  the  109-grid-point  problem  by  using  a 
large-variety  of  high-performance,  expensive  components. 
However,  the  many  different  types  of  components  required 
and  their  diverse  interconnections  would  not  allow  any 


55 


CELLUAR- ARRAY -COMPUTER  PARAMETERS 


tion  of  state  of  the 


regular  construction  pattern,  so  that  the  complexity 
would  be  overwhelmingly  great.  Such  a machine  appears 
impractical  for  the  near  future  (the  next  10  to  15  years). 

(4)  A cellular-array  computer,  composed  of  a single  inexpensive 

general-purpose  minicomputer  driving  a very  large, 

regular  array  of  identical  cells  of  a single  type,  each 

cell  of  which  represents  one  or  more  grid  points,  does 

g 

appear  to  be  feasible  for  the  10  -grid-point  problem 
by  the  early  1980s. 

(5)  A cellular  array  of  more  modest  size  might  be  very 
effective  for  calculating  some  specialized  problems 
in  the  transition  to  turbulence,  plasma  instability, 
etc.  hence,  some  smaller  array  may  be  worthwhile  prior 
to  the  1980s,  when  ISF  technological  advancements 

will  probably  allow  the  construction  of  the  large 
array  for  a cost  that  is  not  totally  unreasonable. 

(6)  A very  modest  investigation  of  the  cellular-array 
computer  appears  to  be  justified,  since  many 
interesting  questions  must  be  answered  before  such 
a machine  could  be  constructed. 

(7)  Extensive  simulation  of  the  array  machine  should 
precede  any  commitment  to  hardware.  This  could  be 
accomplished  by  a continuing  research  program  at 
modest  cost. 

(8)  The  potential  benefits  of  a powerful  computer  of  the 
type  proposed  here  seem  to  be  very  great.  Some 
further  investigation  is  certainly  indicated. 
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A MODEL  TO  STUDY  DOWNSTREAM  BOUNDARY  CONDITIONS 
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Appendix  C 


A MODEL  TO  STUDY  DOWNSTREAM  BOUNDARY  CONDITIONS 


In  the  main  text  a number  of  arguments  were  given  to  suggest  that 
downstream  boundary  conditions  cause  little  effect  at  a few-boundary- 
layer  thickness  upstream.  Here  we  give  a simple  model  in  which  we  can 

verify  the  statements. 

We  consider  Couette  flow  between  a stationary  plate  at  y = 0 and  a 
plate  moving  with  velocity  U in  the  x-direction  at  position  y = 6.  The 
laminar  solution  of  the  Navier-Stokes  equations  is 


v = u (y),  0,  0 , where  u (y)  = Uy/6 

~ o ° 


(Thus  we  are  approximating  the  Blasius-type  profile  by  a straight  line.) 

Further,  we  neglect  viscosity  and  restrict  ourselves  to  flow  linearized 

around  the  laminar  solution.  We  wish  to  see  how  a disturbance  at  a point 

affects  the  flow  upstream  and  downstream  from  that  point.  Writing 

v = (u  + u , v , 0)  we  have  the  equations 
~ oil 

^i 

at  1 uo  jx  i ay  ay 


avx  du 

+ u — 

at  o ax 


-ap/ay 


ax  ay 


0 


Introducing  a stream  function  <p  by 
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acp 

v = - — * 
1 5x 


u 


1 


_ 5cg 

ay 
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we  readily  find  that 


a_ 

at 


+ u 


o ax 


) (§•£)■• 


(Note  that  ^cp/Sx2  + a2co/ay2  is  just  the  vorticity  of  the  perturbed  flow.) 
Taking  ordinary  Fourier  transforms  with  respect  to  space  t nd  one-sided 
ones  in  time — i.e. , 


and 


:,y  ,(L')  = J 


iajt 


cp(x,y,ai)  =J  e cp(x,y,t)  dt  , 
0 


co(k 


,y.uj)  = j" 


_ ^ 

e cp(x,y,uu)  dx 


— we  see  that 


d2S  . 2 w W(k,y ,0) 

2 " k ® = i(ku  -w) 
dy  o 


(C-l) 


where  W(k,y,0)  is  the  spatial  Fourier  transform  of  the  initial  value  of 
the  vorticity. 

Explicitly , 


W(k,y,0) 


Acp(x,y,0)  dk 


If  we  introduce  the  Green's  function  G(j,y/,k)  by 
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G(y,y/,k)  = - 


sin  ky  sin  k(6  - y)  , 

i i—  V £ y 

k sin  k 6 


sin  ky'  sin  k(6  - y)  , 

1 y a y 

k sin  k 6 


the  solution  of  Eq.  'C-l)  is  given  by: 


cp(k,y  ,uu) 


■/ 


o 


Inverting  the  time  Fourier  transform  gives 


cp(k 


,y.t)  = J 


■iku  (y')t 


G(y,y\k)  W(k,y',0)  dy ' 


and  then 


cp(x,y,t) 


ik[x-u  (y')t] 

G(y,y',k)  W(k,y'0)  dk 


(C-2) 


Now,  let  us  consider  the  effect  of  an  initial  disturbance  at  x = 0, 
y = y . Thus, 


Acp(x,y,0)  = 6 (x)  6(y  - y ) 


Then  W(k,y',0)  = 6 (y  - y ) and  Eq.  (C-2)  becomes 

o 


cp(x,y,t) 


ik[x-u  (y  )t] 

° ° G(y,y  , k)  dk 

o 
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The  integral  can  now  be  evaluated  by  residues.  (To  avoid  unnecessary 
repetition  we  report  merely  the  results  for  y < y .)  Thus, 

(a)  Consider  x-u(y)t<0 
o o 

(This  includes  all  upstream  points  and  downstream 
points  after  a sufficiently  long  time.) 

The  result  is: 


® njr/6[x-u  (y  )t]  , n 
cp(x,y,t)  = e 00  (-1) 

n=l 


nn 


sin 


nn\ 


sin 


nit(6  - y ) 
o 


(C-3) 


Notice  that  the  disturbance  falls  off  exponentially 
with  distance  upstream  (in  a length  ~ 6)  and  decreases 
with  time  exponentially  at  any  fixed  point. 

(b)  x - u (y  ) t > 0. 

o o 

These  are  downstream  points  at  sufficiently  small  times. 


<P(x,y,t> 


-nit/6[x-u  (y  )t] 

n=l 


(-1) 

nn 


nny 

sin  — - sin 


nn(6 


y ) 

o 


(C-4) 


Eq. 


Thus  at  a fixed  point  upstream  the  disturbance  grows  [according  to 
(C-4)],  reaches  a maximum,  and  then  decreases  according  to  Eq.  (C-3) 
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Appendix  D 


BOUNDARY  CONDITIONS 


There  is  a general  class  of  problems,  of  which  the  transition  problem 
is  an  example,  that  require  a careful  examination  of  boundary  conditions. 
We  first  examine  the  question  of  what  must  be  specified  in  order  to  arrive 
at  a unique  solution,  and  we  then  examine  a set  of  boundary  conditions 
for  the  specific  problem  under  consideration.  Lastly,  we  show  that  in 
principle  there  is  a computational  scheme  that  follows  from  this  speci- 
fication. 


The  Navier-Stokes  equations  for  an  incompressible  fluid  governing 

two  flows  v and  v are: 

1 2 


°~1  2 2 

— + v . yv  - -yP  + vy  v 

dt  ~1,2  ~1,2  1,2  vv  ~1 ,2 


7 ' Jl,2  " ° ' 


Defining 


v = v 
~1 


v 


2 


P 


P 

2 


We  find 


dv 

— + v • yv 
at  ~1  ~1 


2 

vo  ' = -yl  + vy  v 
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V • v = 0 


~1  * 7~l  " ~2  ’ V~2  ~1  * V(~  + Z'A)  " Z2  * V~2  = ~1  * V~  + Z * Vy2 


Thus, 


Sv 


at  + ^1  * VZ  + Z * VV  = -vP  + vy  V 


V • v = 0 


The  energy  in  the  difference  flow, 


-Ji*2 


dt 


is  governed  by 


at  +J(V  • Ii  r * X • it  - -fv  • Pv  + v/7  . vv  • V 


which  may  be  rewritten 


_ _ 

dt 


Vv  • v + v 

Ci  ^ 


hxi  8v 


- -/(v  • 

r / v2  2\ 

-|ds(Pn«v-vn«v  — + n v —j 

j \ ~ ~ ~ 2 ~ ~i  2 / 


If  we  suppose  that  = v2  on  the  boundaries,  the  surface  terms  vanish. 
Assume  that  yv2  is  bounded  and  its  maximum  is  gi’  en  by 


M = max(yv  ) 
2 
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Then 


— £ 2M(t)E 
ot 


We  can  then  write 


2 J1  M(t  ' )dt  ' 

0 * E(t ) * E(0)e  ° 

from  which  it  follows  that  if  E(0)  = 0,  E(t)  also  vanishes. 

We  have  therefore  proved  that  specifying  the  velocity  on  the  bound- 
aries is  enough  to  provide  a unique  solution. 

We  next  turn  to  examine  the  specification  of  the  boundary  conditions 
appropriate  to  the  problem  at  hand.  The  coordinate  system  is  shown  in 
Figure  D-l.  The  steady  flow  satisfies  the  exact  equations 


z 


BOUNDARY 

LAYER 


y SPACE 
/ FLOW 
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Sv  Sv 

X z 

T — + “ — = ° 

ox  Sz 


The  usual  physical  picture  of  boundary- layer  theory  involves  the  assump- 
tion that  the  flow  is  more  slowly  varying  in  the  x direction  than  in  the 
z direction.  Thus  we  may  write 


a c a 

3x  Sex 


e « l 


and  find 


dv 

X 

dv 

X 

17  ■ ' — 

SP 

/ ^ 
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dv 
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z Sex  z Sz 

= " Sz  + V 1 
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Sv  Sv 

X z 
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The  divergence  equation  clearly  requires 


v = ev 
z z 


which  yields 


dv  Sv  / S2v  , S~v  \ 

x . — x Sr  I x ix| 

Sex  + Vz  Sz  = " Sex  " V I 2C+e  2 / 

\3(ex)  Sz  / 
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In  order  to  balance  the  viscous  and  inertial  forces 


we  then  must  have 


v = vc 

which  finally  yields  the  standard  boundary-layer  equations,  which  we  re 
write  in  the  original  variables 


dv 

X 

j. 

dv 

X 

dP 

d2v 

v x 

dx  z 

dz 

dx 

dz2 

dv 

z 

dv 

z 

o 

d v 

z 

VX  dx  + 

Z Sz 

~ V 

dz2 

dv 

dv 

x z 

T + — = 0 

ox  dz 


Estimating  the  terms  yields 


Sv  d2v 

x ~ x 

x Sx  x 2 
dz 


or 


v 

X 

V ~ — 

x 2 
z 


V 

z 


This  immediately  suggests  the  Blasius  similarity  solution 
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Ml 


v = Uf(Tl) 
x 


with 


T1 


v = -g01> 
z z 

This  leads  to  the  nonlinear  equation 

ff " + 2f = 0 


and  the  boundary  conditions 

f = 0 , f'=0  , z=C 

f ' = 1 z -*  00 

Since  these  solutions  are  only  approximate,  part  of  th^  flow  field  that 
develops  in  time  in  the  computation  will  result  from  the  equation  trying 
to  predict  the  actual  solutions.  We  look  for  solutions  of  the  form 

-iujt  iPy  ifof(x)dx 
. , e e 

f = f(z)e 


and  we  know  that  the  resulting  equation  is  a form  of  the  Orr- Sommer f eld 
equation  which  is  fourth-order  in  z.  T\vo  boundary  conditions  are  provided 
by  the  vanishing  of  the  perturbed  flow  at  the  bottom  and  top,  and  two  more 
involving  the  derivatives  are  provided  by  the  vanishing  of  other  components 
of  the  velocity.  The  resulting  eigenvalue  equation 
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h(u),f5  ,Qf  ,R)  = 0 


ylelds  the  local  spatial  growth.  In  gerer.l,  «l-««  the  boundary- layer 


thickness 


6 ' k vx/U 


and  the  Reynolds  number 


yields  curves  ot  the  tor.  given  in  Figure  D-2.  The  actual  experiment  ot 


GIVEN  FLOW,  X INCREASING 


FIGURE  D-2  STABILITY  DIAGRAM  FOR  BOUNDARY-LAYER  FLOW 

Klebanof f et  al.  (1961)  corresponds  to  point  A tor  the  upstream  condition. 
To  save  calculations!  time,  the  numerical  calculations  may  he  started  at 
point  B,  which  is  100  cm  downstream,  and  the  calculation  can  be  stopped 

at  point  C. 

- _ nhvsically  reasonable  set  of  boundary  conditions  that 
A summary  of  a physica-ixy 

at  least  is  mathematically  unique  is  then  given  by 

(1)  Global  constraint,  J dS  v . n ! 0. 

(2)  Bottom  plate,  v = 0,  no  slip  condition. 


(3)  Top  of  region,  vx  = U(c°) , free  stream. 

(4)  Spanwise,  the  flow  is  symmetric.  Thus,  vy  = 0 at  the 
sides,  and  dvx/^y  and  dvz/dy  both  vanish.  Although  the 
original  uniqueness  proof  was  carried  out  for  specified 
v,  the  surface  terms  vanish  except  for 

2 r I dv  5v  dv  \ 

Jds  " v r =/dx  dz  (vx  s7  * vy  if  * w) 

which  vanishes  for  the  specification  of  derivatives 
given  above . 

(5)  Upstream.  As  discussed  above,  the  conditions  are 

v = U(z)  + 6v 

X X 


v = 6v 

y y 


v = 6 v 
z z 

where  U(z)  is  the  ttlasius  profile  corresponding  to  a 
position  100  cm  dov/nst*-eam  from  the  entrance,  and  6v 
is  a linear  combination  of  twr - and  three-dimensional 
Tollmien-Schlichting  waves.  These  waves  are  chosen 
to  correspond  to  the  driving  frequency  of  the  ribbon 
in  the  Klebanoff  experiment  and  to  have  the  periodicity 
in  y imposed  by  the  spacers. 

(6)  Downstream.  Here  we  assume  that  the  conditions  imposed 
in  (5)  are  simply  carried  downstream  to  a location  at 
x = L corresponding  to  point  C in  the  stability  diagram. 

We  know  that  the  spatial  growth  is  converted  downstream 
with  a velocity  c that  is  less  then  U^.  Thus,  for  some 
period  of  time,  chosen  so  that  transition  occurs  in  the 
volume  of  the  calculational  domain,  the  downstream 
boundary  conditions  will  not  change. 

I 

We  now  wish  to  show  that,  in  principle,  the  specification  given  above 
yields  a method  of  numerical  calculation.  At  t = 0,  the  velocity  field 
is  specified  everywhere  and  for  all  t,  v is  specified  on  the  boundary. 
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From  the  Navier-Stokes  equations  and  V • v = 0,  it  follows  that 


2 

V P = -V 


(s 


A Greens  function  G satisfying 

2 , 

V G = -6 (x  - x ) 


yields  the  solution 


P = 


(Gn 


(v  * Vv)  ds 

r*j  rv 


Since  P is  unknown  on  the  boundary,  choose  n • VG  = 0 there.  We  are 
guaranteed  then  that  G exists.  From  the  normal  component  of  the  Navier- 
Stokes  equation  on  the  surface 


a 2 

— v • n + v • Vv  • n = -n  • VP  + vV  v • n 


we  see  that  n • VP  is  determined  in  terms  of  known  quantities.  Thus, 


= /°[' 


. 2 a 

P = I G vV  v • n - — v • n - v • Vv  • n 

- ~ st  ~ ~ 


rsj 


ds 


" J G(V  • v 


• Vv)  dt 


and  is  determined  everywhere.  Now,  knowing  P we  use 


Sv 

St 


2 

-VP  + vV  v - v 


Vv 


and  the  knowledge  of  v everywhere  to  advance  v forward  in  time  at  every 
interior  point  of  the  calculational  mesh.  In  a crude  sense  this  represents 
a proof  of  the  existence  of  solutions  to  these  equations. 
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Appendix  E 


THREE-DIMENSIONAL  IN-PLACE  FAST- FOURIER  TRANSFORM  ARRAY 


Preceding  page  blank 
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Appendix  F 


THREE  DIMENSIONAL  IN-PLACE  FAST-FOUR I ER- TRANSFORM  ARRAY 

Assume  that  it  is  desired  that  an  array  of  cePs  is  to  be  permanently 

connected  with  as  few  connections  as  possible  such  that  three-dimensional, 

in-place,  fast  Fourier  transforms  can  be  made  of  variables  stored  within 

the  cells.  That  is,  if  the  elements  of  matrix  A,  a are  stored  in 

i jk 

Cellijk’  thGn  matrix  B w^h  elements  b^  are  to  be  calculated  according 
to 


B = F[A] 


L-l  M-l  N-l 


4 m 


exp 


111  im 
2jf  \ L + M 


Stone  (1970)  has  shown  how  a one-dimensional  FFT  on  an  array  of  cells  can 
be  accomplished  in  place  by  a "perfect  shuffle"  interconnection  between 
cells.  Only  two  inputs  and  two  outputs  per  cell  are  required.  Despain 
(1974)  has  shown  now  FFT  operations  such  as  those  used  within  the  cells 
used  by  Stone  can  be  implemented  without  calculating  or  storing  the 
trigonometric  coefficients  usually  employed  (cor die  methods).  This  makes 
practical  an  array  of  simple  identical  cells  that  can  be  made  to  perform 
an  in-place  transform  of  variables  stored  in  the  cells. 

Three-dimensional  transforms  are  generally  accomplished  as  three 
successive  groups  of  one -dimensional  transforms.  However,  it  is  possible 
to  simultaneously  transform  in  all  three  dimensions  at  the  same  time. 

This  is  of  great  advantage  for  the  three-dimensional  array  of  cells,  since 
the  perfect  shuffle"  connections  are  permanently  connected  in  all  three 

Preceding  page  blank 
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dimensions  and  parallel  data  transfer  (and  calculations  too  if  desired) 
can  proceed  simultaneously  and  at  a much  greater  rate  than  if  separate 
transforms  are  taken  in  succession. 

The  connection  pattc  n is  symmstric  in  all  three  dimensions.  Figures 
E-l  and  E-2  illustrate  these  patterns  for  one  and  two  dimensions  . 


FIGURE  E-1  PERFECT  SHUFFLE  CONNECTIONS 
FOR  1-D  FFT,  N = 16 


FIGURE  E-2  4-BY-4  ARRAY  OF  CELLS  FOR  A 

2-D  FFT  NETWORK 
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